Let’s say you do all that and look at your model output and the estimate is…not significant!
Specify our model \[ \operatorname{Work} = \alpha + \beta_{1}(\operatorname{Coffee}) + \epsilon \] Where, the average effect of coffee, \(beta_1\), is set at 0.1, and \(\epsilon\) is a random error, normally distributed with mean 0, sd of 1. If our model is correctly specified, we should recover this effect.
Let’s look at the residuals of our model along with the spatial lag of residuals (i.e., the residuals of nearest neighbors).
Moran’s I test “measures spatial autocorrelation based on both feature locations and feature values simultaneously.” Basically, it is a measure of how similar each unit’s residual is with some \(k\) set of neighbors.
In math: \[ I = \frac{N}{W} \cdot \frac{\sum_{i=1}^{N} \sum_{j=1}^{N} w_{ij} (x_i - \bar{x})(x_j - \bar{x})}{\sum_{i=1}^{N} (x_i - \bar{x})^2} \] Where, \(N\) is the sample size for unit \(i\) across area \(j\), \(x\) is the variable of interest, \(W\) is the weight created through row-wise standardization so that the sum of all of unit \(i\)’s neighbor’s weights is equal to 1. The ArcGIS documentation on this is very good.
Let’s start with a test that assesses the correlation with a given unit’s 15 nearest neighbors. Moran’s I is based on a basic hypothesis testing framework, where:
A positive Moran’s I suggests positive clustering, i.e., units tend to have similar values, while a negative value suggests dispersion different from random.
Global Moran I for regression residuals
data:
model: lm(formula = Work ~ Coffee + Age + Income, data = data)
weights: listw
Moran I statistic standard deviate = 13.552, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
0.422427205 -0.010046388 0.001018314
When we talk about p-values, what we are talking about is the probability of obtaining a Moran’s I at least as extreme as the one observed assuming the null is true for our sample. The null implies random distribution of residuals, or the middle box below.
The spatial error model is essentially adding a weighted residual to account for the fact that there are relationships not controlled for in your original specification.
\[ {\mathbf y} = {\mathbf X}{\mathbf \beta} + {\mathbf u}, \qquad {\mathbf u} = \rho_{\mathrm{Err}} {\mathbf W} {\mathbf u} + {\mathbf \varepsilon} \]
Put simply, sometimes the omitted variable biasing your results is right there next to you, and your neighbor, and your neighbor’s neighbor.
This is a map of depression in Seattle, Washington. We want to understand the relationship between neighborhood characteristics and depression prevalence.
We can map the residuals of a basic linear model to get a better idea of correlation
Global Moran I for regression residuals
data:
model: lm(formula = DEP_CrudePrev ~ unempr + pmob + pcol + ppov +
pnhblk + phisp + log(tpop), data = sea.tracts)
weights: seaw
Moran I statistic standard deviate = 6.6676, p-value = 1.3e-11
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
0.309531240 -0.026084108 0.002533614
Thank you!